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Abstract — Emerging sonography techniques often require in- 
creasing the number of transducer elements involved in the 
imaging process. Consequently, larger amounts of data must be 
acquired and processed. The significant growth in the amounts 
of data affects both machinery size and power consumption. 
Within the classical sampling framework, state of the art systems 
reduce processing rates by exploiting the bandpass bandwidth of 
the detected signals. It has been recently shown, that a much 
more significant sample-rate reduction may be obtained, by 
treating ultrasound signals within the Finite Rate of Innovation 
framework. These ideas follow the spirit of Xampling, which 
combines classic methods from sampling theory with recent 
developments in Compressed Sensing. Applying such low-rate 
sampling schemes to individual transducer elements, which detect 
energy reflected from biological tissues, is limited by the noisy 
nature of the signals. This often results in erroneous parameter 
extraction, bringing forward the need to enhance the SNR of the 
low-rate samples. In our work, we achieve SNR enhancement, 
by beamforming the sub-Nyquist samples obtained from multiple 
elements. We refer to this process as "compressed beamforming". 
Applying it to cardiac ultrasound data, we successfully image 
macroscopic perturbations, while achieving a nearly eight-fold 
reduction in sample-rate, compared to standard techniques. 

Index Terms — Array Processing, Beamforming, Compressed 
Sensing (CS), Finite Rate of Innovation (FRI), Ultrasound, 
Xampling 

I. Introduction 

Diagnostic sonography allows visualization of body tissues, 
by radiating them with acoustic energy pulses, which are 
transmitted from an array of transducer elements. The image 
typically comprises multiple scanlines, each constructed by 
integrating data collected by the transducers, following the 
transmission of an energy pulse along a narrow beam. As 
the pulse propagates, echoes are scattered by density and 
propagation- velocity perturbations in the tissue and de- 
tected by the transducer elements. Averaging the detected 
signals, after their alignment with appropriate time- varying 
delays, allows localization of the scattering structures, while 
improving the Signal to Noise Ratio (SNR) |2|. The latter 
process is referred to as beamforming. Performed digitally, 
beamforming requires that the analog signals, detected by the 
transducers, first be sampled. Confined to classic Nyquist- 
Shannon sampling theorem |3|, the sampling rate must be at 
least twice the bandwidth, in order to avoid aliasing. 

As imaging techniques develop, the amount of elements 
involved in each imaging cycle typically increases. Conse- 
quently, the rates of data which need to be transmitted from the 
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system front-end, and then processed by the beamformer, grow 
significantly. The growth in transmission and processing rates 
inevitably effects both machinery size and power consumption. 
Consequently, in recent years there has been growing interest 
in reducing the amounts of data as close as possible to the 
system front-end. In fact, such reduction is already possible 
within the classical sampling framework: state of the art 
devices digitally downsample the data at the front-end, by 
exploiting the fact that the signal is modulated onto a carrier, 
so that the spectrum essentially occupies only a portion of 
its entire base-band bandwidth. The preliminary sample rate 
remains unchanged, since the demodulation is performed in 
the digital domain. Nevertheless, a key to significant data 
compression lies beyond the classical sampling framework. 

Indeed, the emerging Compressive Sensing (CS) frame- 
work |4|, |5| states, that sparse signals may be accurately 
reconstructed from a surprisingly small amount of coefficients. 
Complementary ideas rise from the Finite Rate of Innovation 
(FRI) framework |6|, in which the signal is assumed to 
have a finite number of degrees of freedom per unit time. 
Many classes of FRI signals can be recovered from samples 
taken at the rate of innovation fT|. For a detailed review of 
previously proposed FRI methods, the reader is referred to 0. 
Combining the latter notions with classical sampling methods, 
the developing Xampling framework f9l, fTOl, flT| involves 
methods for fully capturing the information carried by an 
analog signal, by sampling it far below the Nyquist-rate. 

Following the spirit of Xampling, Tur et. al. proposed 
in [121, that ultrasound signals be described within the FRI 
framework. Explicitly, they assume that these signals, formed 
by scattering of a transmitted pulse from multiple reflectors, 
may be modeled by a relatively small number of pulses, all 
replicas of some known pulse shape. Denoting the number of 
reflected pulses by L, and the signal's finite temporal support 
by [0,T), the detected signal is completely defined by 2L de- 
grees of freedom, corresponding to the replicas' unknown time 
delays and amplitudes. Based on |6|, the authors formulate the 
relationship between the signal's Fourier series coefficients, 
calculated with respect to [0, T), and its unknown parameters, 
in the form of a spectral analysis problem. The latter may be 
solved using existing techniques, given a subset of Fourier 
series coefficients, with a minimal cardinality of 2L. The 
sampling scheme is thus reduced to the problem of extracting a 
small subset of the detected signal's frequency samples. Two 
robust schemes are derived in lITZl . llT3l , extracting such a 
set of coefficients from samples of the signal, taken at sub- 
Nyquist rates. The system presented in llT2l employs a single 
processing channel, in which the analog signal is filtered 
by an appropriate sampling kernel and then sampled with 
a standard low-rate analog to digital converter (ADC). The 
method of IIT3l employs multiple processing channels, each 
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comprising a modulator and an integrator. These approaches 
were shown to be more robust than previous FRI techniques 
and also allow for arbitrary pulse shapes. 

The initial motivation for our work stems from the need to 
translate the ultrasound Xampling scheme proposed in |12], 
into one which achieves the final goal of reconstructing a two- 
dimensional ultrasound image, by integrating data sampled at 
multiple transducer elements. In conventional ultrasound imag- 
ing, such integration is achieved by the beamforming process. 
The question is how may we implement beamforming, using 
samples of the detected signals taken at sub-Nyquist rates. 

A straightforward approach is to replace the Nyquist-rate 
sampling mechanism, utilized in each receiver element, by an 
FRI Xampling scheme. Having estimated the parametric rep- 
resentation of the signal detected in each individual element, 
we could reconstruct it digitally. The reconstructed signals 
can then be further processed via beamforming. However, 
the nature of ultrasound signals reflected from real tissues, 
makes such an approach impractical. This is mainly due to 
the detected signals' poor SNR, which results in erroneous 
parameter extraction by the Xampling scheme, applied to each 
element independently. 

Our approach is to generalize the FRI Xampling scheme 
proposed in |13|, such that it integrates beamforming into 
the low-rate sampling process. The result is equivalent to that 
obtained by Xampling the beamformed signal, which exhibits 
significantly better SNR. Furthermore, beamforming practi- 
cally implies that the array of receivers is dynamically focused 
along a single scanline. Consequently, the resulting signal 
depicts reflections originating in the intersection of the radiated 
medium with a vary narrow beam. Such a signal better suits 
the FRI model proposed in |T2], which assumes the reflections 
to be caused by isolated, point-like scatterers. We refer to our 
scheme by the term compressed beamforming, as it transforms 
the beamforming operator into the compressed domain [14], 
idgl. Applied to real cardiac ultrasound data obtained from a 
GE breadboard ultrasonic scanner, our approach successfully 
images macroscopic perturbations in the tissue while achieving 
a nearly eight-fold reduction in sampling rate, compared to 
standard imaging techniques. 

The paper is organized as follows: in Section [III we sum- 
marize the general principles of beamforming in ultrasound 
imaging. In Section [nil we outline the FRI model and its 
contribution to sample rate reduction in the ultrasound con- 
text. We motivate compressed beamforming in Section HVl 
considering the nature of ultrasound signals reflected from 
biological tissues. Beamforming and FRI Xampling are com- 
bined in Section |Vl where we propose that the signal obtained 
by beamforming may be treated within the FRI framework. 
Following this observation, we derive our first compressed 
beamforming scheme, which operates on low-rate samples 
taken at the individual receivers. This approach is then further 
simplified in Section [Vll In Section IVIII we focus on image 
reconstruction from the parametric representation obtained by 
either Xampling scheme. In this context, we generalize the 
signal model proposed in [12], allowing additional unknown 
phase shifts of the detected pulses. We then discuss an alterna- 
tive recovery approach, based on CS. Simulations comparing 
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Fig. 1. Imaging setup: M receivers are aligned along the x axis. The origin 
is set at the position of the reference receiver, denoted mo- Sm denotes 
the distance measured from the reference receiver to the mth receiver. The 
imaging cycle begins when an acoustic pulse is transmitted at direction 6. 
Echoes are then reflected from perturbations in the radiated medium. 

the performance of several recovery methods are provided in 
Section rvllll Finally, experimental results obtained for cardiac 
ultrasound data are presented in Section [1x1 

II. Beamforming in Ultrasound Imaging 

In this section, we describe a typical B-mode imaging cycle, 
focusing on the beamforming process, carried out during the 
reception phase. The latter constitutes a significant block in 
ultrasound imaging, and plays a major role in our proposed 
FRI Xampling scheme. 

Consider the array depicted in Fig. [T] comprising M trans- 
ducer elements, aligned along the x axis. Denote by 5m the 
distance from the mth element to the reference receiver mo, 
used as the origin, namely 5 mo = 0. The imaging cycle begins 
when, at time t = 0, the array transmits acoustic energy into 
the tissue. Subsequently, the elements detect echoes, which 
originate in density and propagation- velocity perturbations, 
characterizing the radiated medium. Denote by (pm {t) the 
signal detected by the mth receiver. The acoustic reciprocity 
theorem |16| suggests, that we may use the signals detected 
by multiple transducer elements, in order to probe arbitrary 
coordinates for reflected energy. Namely, by combining the 
detected signals with appropriate time delays, echoes scattered 
from a chosen coordinate will undergo constructive interfer- 
ence, whereas those originating off this coordinate will be 
attenuated, due to destructive interference. 

In practice, the array cannot effectively radiate the entire 
medium simultaneously. Instead, a pulse of energy is con- 
ducted along a relatively narrow beam, whose central axis 
forms an angle with the z axis. Focusing the energy pulse 
along such a beam is achieved by applying appropriate time 
delays to modulated acoustic pulses, transmitted from multiple 
array elements. Rather than arbitrarily probing the radiated 
tissue, we are now forced to adjust the probed coordinate in 
time, in coordination with the propagation of the transmitted 
energy. This practically implies that, combining the detected 
signals with appropriate time- varying delays, we may obtain a 
signal, which depicts the intensity of the energy reflected from 
each point along the central transmission axis. Throughout 
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the rest of this section, we derive an expHcit expression for 
creating this beamformed signal. 

Assume that the energy pulse, transmitted at t = 0, 
propagates at velocity c in the direction 0. At time t > 0, 
the pulse crosses the coordinate (x^z) = {ct sin ^ ct cos 0) . 
Consider a potential reflection, originating in this coordinate, 
and arriving at the mth element. The distance traveled by such 
a reflection is: 

dm{t; 9) = \J {ct cos + {Sm - ct sin Of . (1) 

The time in which the reflection crosses this distance is 
dm {t; 0) /c, so that it reaches the receiver element at time 

f^(i;^)=t+^!!iM. (2) 

c 

It is readily seen that Tmo {t; 0) = 2t. Hence, in order to 
align the reflection detected in the mth receiver with the one 
detected in the reference receiver, we need to apply a delay 
to (pm{t), such that the resulting signal, (prn{t'-,0), satisfies 
(fm (2t; 0) = (fm {fm [t] 6)). Denoting (t; 9) = fm {t/2; 0), 
and using ([T]), we obtain the following distorted signal for 
t > 0: 

rm {t; ^) = ^ + ^yt^-4jmtsm0^4jl^ , 

with = ^m/c- The aligned signals may now be averaged, 
resulting in the beamformed signal 

1 ^ 

^it;0) = -Y.^rn{t;0), (4) 

m=l 

which exhibits enhanced SNR compared to {(pm {t;0)}^^-^. 
Furthermore, by its construction, ^ {t;6) represents, for every 
t > 0, the intensity which was measured when focusing the 
array to p{t) = {ct/2sm0^ct/2cos6). Therefore, it may 
eventually be translated into an intensity pattern, plotted along 
the corresponding ray. 

Although defined over continuous time, ultrasound systems 
perform the process formulated in (O-© in the digital domain, 
requiring that the analog signals (pm {t) first be sampled. 
Confined to the classic Nyquist- Shannon sampling theorem, 
these systems sample the signals at twice their baseband 
bandwidth, in order to avoid aliasing. The detected signals 
typically occupy only a portion of their baseband bandwidth. 
Exploiting this fact, some state of the art systems manage 
to reduce the amount of samples transmitted from the front- 
end, by down- sampling the data, after demodulation and low- 
pass filtering. However, since such operations are carried out 
digitally, the preliminary sampling-rate remains unchanged. 

To conclude this section, we evaluate the nominal number of 
samples needed to be taken from each active receiver element 
in order to obtain a single scanline using standard imaging 
techniques. Consider an ultrasound system which images to a 
nominal depth of r = 16cm. The velocity at which the pulse 
propagates, c, varies between 1446m/sec (fat) to 1566m/sec 
(spleen) ifTTl . An average value of 1540m/sec is assumed by 
scanners for processing purposes, such that the duration of the 



detected signal is T = 2r/c ^ 210/isec. The signal's baseband 
bandwidth requires a nominal sampling rate of fs = 16Mhz, 
resulting in an overall number of Tfs = 3360 real-valued 
samples. Assuming that the signal's passband bandwidth is 
only 4MHz, the data sampled at Nyquist-rate may be finally 
down-sampled to 1680 real- valued samples. These samples, 
taken from all active receivers, are now processed, according 
to (O-©, in order to construct the beamformed signal. Since 
standard imaging devices carry out beamforming by applying 
delay and sum operations to the sampled data, the amount of 
operations required for generating a single scanline is directly 
related to the sample rate. 

Regardless of our computational power, physical constraints 
imply that the time required for constructing a single scanline 
is at least T. This takes into account the round-trip time re- 
quired for the transmitted pulse to penetrate the entire imaging 
depth, and for the resulting echoes to cross a similar distance 
back to the array. Nevertheless, sufficient computational power 
may allow construction of several scanlines, within that same 
time interval, increasing the overall imaging rate. By using 
compressed beamforming, we aim at capturing significant 
information in the imaging plane, while reducing the sampling 
rate and consequently the processing rate. This, in turn, may 
improve the existing trade-off between imaging rates and both 
machinery size and power consumption. 

III. Sample Rate Reduction Using the FRI Model 

In a pioneer attempt to implement Xampling methodology 
in the context of ultrasound imaging, [12J suggests that the 
signal detected in each receiver element may be sampled at a 
rate far below Nyquist, by modeling it as an FRI signal. The 
authors propose that cpm {t), detected in the mth element, be 
regarded as sum of a relatively small number of pulses, all 
replicas of some known pulse shape. Explicitly: 

L 

(t) = ^ ai^rnh {t - ti^rn)- (5) 

Here L is the number of scattering elements, distributed 
throughout the sector radiated by the transmitted pulse, U^rn 
denotes the time in which the reflection from the /th element 
arrived at the mth receiver, and ai^rn denotes the reflection's 
amplitude, as detected by the mth receiver. Finally, h (t) 
denotes the known pulse shape, regarded, in our work, by the 
term two-way pulse. The signal in © is completely defined 
by 2L real-valued parameters, {ti^rnjCii,m}f^i- 

Sampling FRI signals was first treated by Vetterli et. al. 
|6|. Their approach involves projecting the FRI signal, char- 
acterized by 2L degrees of freedom per unit time, onto a 2L- 
dimensional subspace, corresponding to a subset of its Fourier 
series coefficients. Having extracted 2L frequency samples of 
the signal, spectral analysis techniques (e.g. annihilating fil- 
ter ifTSl . matrix pencil 1 19|) may be applied, in order to extract 
the unknown signal parameters. Applying this solution to the 
problem formulated in (O, flT] formalizes the relationship 
between the ultrasound signal's Fourier series coefficients to 
its unknown parameters, as a spectral analysis problem. 
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Let T be the duration of (^rn {t)- We can then expand (pm (t) 
in a Fourier series, with coefficients 

1 

^ Jo 
1 ^ 

^cii.mh(t-ti^m)e~'^^^dt (5) 
1=1 

where H {uj) denotes the Continuous Time Fourier Transform 
(CTFT) of h (t). Consider the sequence {^j,m}^^? compris- 
ing Km integers, and define the length-K^ vector with 
jth element (j)rn [^j,m]- Then ^ may be written in matrix 
form: 

^rn = ^HmVmam, (7) 

where Hm is a diagonal matrix with diagonal elements 
H {^kj^rn), Vm contains e"*^^^^^ "^ as its {jjl)th element, 
and am is the length L vector, with elements ai^rn- Choosing 
kj^rn such that H {^kj^rn) 7^ 0,we can express © as: 

Ym = V^am, (8) 

where ym = TH^^^^m- If the values kj^rn are a sequence 
of consecutive indices, then Vm takes on a Vandermonde 
form, and has full column rank | IH as long as Km > L and 
the time-delays are distinct, i.e., ti^m 7^ ^j,m, for all i ^ j. 
The formulation derived in ([5]) is a standard spectral analysis 
problem. As long as Km > 2L, it may be solved for the 
unknown parameters {ti,m^CLi,m}f=i, using methods such as 
annihilating filter 1 1 8 1 or matrix pencil L19J . 

Having obtained (|7]), the sampling scheme reduces to the 
problem of extracting Km frequency samples of ipm {t), where 
Km > 2L. A single-channel Xampling scheme, such as the 
one derived in |[T2ll . allows robust estimation of such coeffi- 
cients from point-wise samples of the signal, after filtering it 
with an appropriate kernel. The estimation is performed by 
applying a linear transformation to p complex- valued samples 
(equivalently, 2p real-valued samples) of the filtered signal, 
requiring that p > Km- In this context, [12] introduces the 
Sum of Sines kernel, which satisfies the necessary constraints, 
and is additionally characterized by a finite temporal support. 
Combining the requirements that Km >2L and p > Km, the 
Xampling scheme proposed in lfT2l allows reconstruction of 
the signal detected in each receiver element from a minimal 
number of 41/ real- valued samples. Considering the nominal 
figures derived in the previous section for standard beamform- 
ing, we conclude that, as long as 4L <C 1680, such a Xampling 
method may indeed achieve a substantial rate reduction. 

IV. Why Compressed Beamforming? 

Applied to a single receiver element, the Xampling scheme 
proposed in llT2l achieves good signal reconstruction for an 
actual ultrasound signal, reflected from a setup of phantom 
targets. In principle, we could apply this approach to each 
receiver element individually, resulting in a^arametric rep- 
resentation for each of the signals {(pm {t)}^^^. Being able 



to digitally reconstruct the detected signals, we could then 
proceed with the standard beamforming process, outlined in 
Section [III aimed at constructing the corresponding scanline. 
Computational effort would have been reduced, by limiting 
the beamforming process to the support of the estimated 
pulses. In fact, we could possibly bypass the beamforming 
stage, by deriving a geometric model which maps the set 
of delays, {ti^m}m=v associated with the Ith reflector, to 
its two-dimensional position pi = {xi^zi). However, apply- 
ing the proposed FRI Xampling scheme to signals reflected 
from biological tissues, we face two fundamental obstacles: 
low SNR and proper interpretation of the estimated signal 
parameters, considering the profile of the transmitted beam. 
These two difficulties may be better understood by examining 
Fig. [21 which depicts traces acquired for cardiac images of a 
healthy consenting volunteer using a GE breadboard ultrasonic 
scanner. 

In the left plot (a), are signals detected by 32 of 64 active 
array elements, following the transmission of a single pulse. 
The pulse was conducted along a narrow beam, forming an 
arbitrary angle with the z axis. The right plot (b) depicts 
the signal obtained by applying beamforming to the detected 
signals, as outlined in Section III Examining the individual 
traces, one notices the appearance of strong pulses, possibly 
overlapping, characterized by a typical shape, as proposed in 
(0). Let us assume that we could indeed extract the delays 
and amplitudes of these pulses, by applying the proposed FRI 
Xampling scheme to each element. We suggested that beam- 
forming could be bypassed, by deriving a geometric model 
for estimating the two-dimensional position of a scattering 
element, based on the delays of pulses associated with it, 
yet estimated in different receivers. In order to apply such 
a model, we must first be able to match corresponding pulses 
across the detected signals. However, referring to the practical 
case depicted in (a), we notice that such a task is not at all 
trivial - the individual signals depict reflections, originating 
from the entire sector, radiated by the transmitted pulse. These 
reflections may, therefore, vary significantly across traces. In 
fact, some pulses, visible in several traces, are not at all 
apparent in other traces. In contrast, the beamformed signal, 
by its construction, depicts intensity of reflections originating 
from along the central transmission axis, while attenuating 
reflections originating off this axis. 

Attempting to apply FRI Xampling to each receiver element 
individually, we encounter an even more fundamental obsta- 
cle, at the earlier stage of extracting the signal's parametric 
representation from its low-rate samples. The individual traces 
contain high levels of noise. The noisy components, especially 
noticeable in traces 54 and 64, rise mainly from constructive 
and destructive interference of acoustic waves, reflected by 
dense, sub-wavelength scatterers in the tissue. The latter are 
typically manifested as granular texture in the ultrasound im- 
age, called speckle, after a similar effect in laser optics L2J . The 
noisy components inherently induce erroneous results, when 
attempting to sample and reconstruct the FRI components 
using the Xampling approach. In extreme scenarios, where 
the noise masks the FRI component, the extracted parameters 
will be meaningless, such that any attempt to cope with errors 
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in the parametric domain will turn out useless. 

The motivation to our approach rises from the observa- 
tion, that we may resolve the aforementioned obstacles by 
Xampling the beamformed signal, ^{t;0), rather than the 
individual signals (pm {t)- Whereas beamforming is a funda- 
mental process in ultrasound imaging since its early days, our 
innovation regards its integration into the Xampling process. 
We derive our compressed beamforming approach, beginning 
with conceptual Xampling of the beamformed signal, using 
the scheme proposed in (13] . We then show that an equivalent 
result may be obtained from low-rate samples of the individual 
signals (pm {t)- 

A necessary condition for implementing our approach is 
that ^{t]0), generated from {^m{t)}^^i satisfying (O, is 
also FRI of similar form. Examining Fig. |2] we notice that 
^ {t;0) exhibits a structure similar to that of the individual 
signals, comprising strong pulses of typical shape, which may 
overlap. In this case, there are several obvious advantages in 
Xampling ^ {t;0). First, since {(pm (^)}^=i are averaged in 
^ {t;6) (after appropriate distortion, derived from the acoustic 
reciprocity theorem) it naturally exhibits enhanced SNR with 
respect to the individual signals. The attenuation of noise in 
the beamformed signal, compared to the individual signals, is 
apparent in Fig. O especially in the interval 50mm — 80mm. 
Second, ^ {t;6) is directly related to an individual scanline. 
This means that we are no longer bothered with the ambiguous 
problem of matching pulses across signals detected in different 
elements. Finally, recall that the signal model derived in (O 
assumes isolated point-reflectors. Such a model is better justi- 
fied with respect to ^ {t; 0) since, by narrowing the effective 
width of the imaging beam, we may indeed approximate its 
intersection with reflecting structures to be point-like. This 
effect is noticeable in Fig. [2] where some pulses, visible in 
individual traces, appear attenuated in the beamformed signal. 
Such pulses correspond to reflectors located off the central 
axis of the transmission beam. 

In the next section, we focus on justifying the assumption 
that ^ {t]0) may be treated within the FRI framework. An 
additional challenge, implied in Section HIl regards the fact 



that ^ {t;6) does not exist in the analog domain - standard 
ultrasound devices generate it digitally, from samples of the 
signals detected in multiple receiver elements, taken at the 
Nyquist-rate. Our goal is, therefore, to derive a scheme, which 
manages to estimate the necessary samples of ^^;6), from 
low-rate samples of filtered versions of {(pm {t)}m=i- 

V. Compressed Beamforming 

Our approach is based on the assumption that the FRI 
scheme, outlined in Section |III1 may be applied to the beam- 
formed signal ^{t;0), constructed according to (O-©. The 
latter exhibits much better SNR than signals detected in 
individual receiver elements. Additionally, it depicts reflections 
originating from a sector much narrower than the one radiated 
by the transmission beam. Its translation into a single scanline 
is therefore straightforward. In Section IV-AI we prove that if 
the signals (pm {t) obey the FRI model (0, then ^ {t;6) is 
approximately of the form: 

L 

^{t',e) = Y,hih{t-ti), (9) 

where ti denotes the time in which the reflection from the lih 
element arrived at the reference receiver, indexed mo. ^ {t;0) 
may thus be sampled using the Xampling schemes derived 
in [if2ll . ifTSl . In practice, we cannot sample <l> {t; 0) directly, 
since it does not exist in the analog domain. In Second IV-BI 
we show how the desired low-rate samples of ^ [t] 6) can be 
determined from samples of (prn {t)- 

A. FRI Modeling of the Beamformed Signal 

Throughout this section we apply three reasonable assump- 
tions. First, we assume that 27^ < ti. Practically, such 
a constraint may be forced by appropriate apodization, as 
often performed in ultrasound imaging. Namely, (prn (f) is 
combined in $ (t]0) only for t > 27^. As an example, for 
the breadboard ultrasonic scanner used in our experiments, the 
array comprised 64 receiver elements, distanced 0.29mm apart. 
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The proposed apodization implies that the receivers located 
farthest from the origin are combined in the beamformed 
signal for imaging depth greater than 9.1mm. Second, we 
assume the two-way pulse, h{t), to be compactly supported on 
the interval [0, A). Finally, we assume that A <C t^. The last 
assumption may also be forced by appropriate apodization. As 
an example, the nominal duration of the pulse acquired by the 
breadboard ultrasonic scanner used in our experiments was 
4/isec. In this case, echoes scattered from depth greater than 
3.1cm already satisfy ti > 10 A. 

Suppose that (fm {t) can be written as in (|5]). Applying the 
beamforming distortion ([3]), we get 



where 



E 

1=1 



h (t„ {t; 9) - U 



(10) 



The resulting signal comprises L pulses, which are distorted 
versions of the two-way pulse h{t). Suppose that some of 
the pulses originated in reflectors located off the central beam 
axis. Beamforming implies that, once averaging the distorted 
signals according to d?]), such pulses will be attenuated due 
to destructive interference. Being interested in the structure of 
the beamformed signal ^{t;0), we are therefore concerned 
only with pulses which originated in reflectors located along 
the central beam. For convenience, we assume that all pulses 
in ([Tot satisfy this property (pulses which do not satisfy it, 
will vanish in ^ {t;6)). We may thus use Tm {t;0), defined in 
©, in order to express U^rn in terms of U. Substituting t = ti 
into Tm {t; 0), we get ti^m = Tm {U; 0), so that (fTOl) becomes 



1=1 



(11) 



where we defined hi^rn {t'l 0) = h {rm {t; 0) — {ti;0)). 

Applying our second assumption, the support of hi^rn {t; 0) 
is defined by the requirement that 



{)<T^(t-e)-T^(ti-e)</\. 



(12) 



Using (O and ©, it is readily seen that hi^rn {t; 0) is 
supported on [U^ti + A'), where 



A' = 2A 



^tj-^lmti sin i9 + 47^ + A 



- ^Imti sin i9 + 472^ + 2 A + - 27^ sin 6 ' 

(13) 



Further applying our assumption that 27^ < tu we obtain 
A' < 2A. 

We have thus proven that hi^rn {t'-, 0) = for t ^ 
[ti,ti^2A). Next, let us write any t in [ti,ti^2A) as 
t = ti^r], where < 7^ < 2A. Then 

hm {t; 0) = h {Tm {tl -^mO)- Tm (tl'.O)) . (14) 

We now rely on our assumption that A <C t/. Since r] < 2 A, 
we also have r] <^ ti. The argument of h{-) in (O may 
therefore be approximated, to first order, as 

Tm {tl + 6>) - Tm {U; 0) = am,l {0) T] ^ O {f]^) , (15) 



^m,l (^) = 2 1 



tl - 2-fm sin 



^if-^lmti sin 19 + 472^ 



(16) 



Up until this point, we assumed that 27^ < ti. Further 
assuming that 7^ <C crm,i {0) 1. Replacing r] by 
T] = t — tl, ([T4l) may therefore be written as 



hm{t;0)^h{t-ti) te [ti.ti ^2A). 



(17) 



Combining (fTTl) with the fact that h{t — ti) is zero outside 
[ti^ti + 2A), (fTTl) may be approximated as 

L 

^m{t;0)^Y.ai^mh{t-ti), (18) 

1=1 

Averaging the signals {(pm{t; 0)}^^-^ according to Q, we get: 

' 1 ^ ^ 



h{t-ti) = Yhih{t-ti), 



=1 



(19) 

which is indeed the FRI form (|9]). Additionally, assuming 
that the support of (^rn (f) is contained in [0,T), we show 
in the Appendix that there exists Tb iO) < T, such that the 
support of <l> (t; 0) is contained in [0, Tb iO)) and, additionally, 

Tm{TB{e)',e)<T. 

As 7m grows towards tu crm,i {0) decreases, resulting in 
a larger distortion of the Ith pulse. Consequently, the ap- 
proximation of (f)rn{t]9) as a sum of shifted replicas of the 
two-way pulse becomes less accurate. The Xampling schemes 
used by lfT2ll . ifTSl rely on the projection of the detected 
signal onto a subspace of its Fourier series coefficients. We 
therefore examine the dependency of the projection error 
on the distortion parameters, 7^, ti and 0. In Fig. [3l we 
show projection errors calculated numerically, for a signal 
comprising a single pulse of duration A = 2/isec. The pulse 
was simulated by modulating a Gaussian envelope with carrier 
frequency 3MHz. It was then shifted by multiple time delays, 
tu where < ti < T, and T = 210/isec, corresponding 
to an imaging depth of 16cm. For each delay, we generated 
the signals (pm {t), assuming that the reflector is positioned 
along the z axis (6 = 0), and that the receiver elements are 
distributed 0.29mm apart, along the x axis. We chose M = 63, 
such that the center (reference) receiver was indexed mo =32. 
The beamforming distortion was then applied to the simulated 
signals, based on ([3]). Finally, the distorted signals were 
projected onto a subset of iiT = 121 consecutive Fourier series 
coefficients, taken within the essential spectrum of the two- 
way pulse. The coefficients extracted from the mth distorted 
signal were arranged into the length K vector, As implied 
by ([3]), no distortion is applied to the signal detected at the 
reference receiver. We therefore evaluate the projection error 
by calculating the SNR defined as 201ogio J^^^^^f-J^- 

The traces obtained for several values of 1 < m < 32 
are depicted in the figure. As ti grows, am,i {0) approaches 
1, and the approximation ([TSl) becomes more valid. As a 
result, the projection error decreases. For receivers located 
near the origin, such that Sm <C 10mm, the error decreases 



Fig. 3. Projection error caused by beamforming distortion with ^ = vs. 
pulse delay, ti, for several receiver elements. The elements are distributed 
0.29mm apart, such that 5i = 8.99mm (element farthest from array center) 
and 5^1 = 0.29mm. Zero error is obtained for the center element, ^32, since 
no distortion is required in this case. 



very quickly. For instance, examining J31 = 0.29mm, the 
SNR grows above 25dB for a reflection originating at distance 
greater than 1/50 of the imaging depth. The SNR improves 
more moderately for receivers located farther away from the 
origin. Nevertheless, considering the receiver located farthest 
away from the origin, ^1 = 8.99mm, the SNR grows above 
lOdB for a reflection originating at distance greater than 1/5 
of the imaging depth. 

Concluding this section, our empirical results indeed justify 
the approximation proposed in Q, where appropriate apodiza- 
tion may further improve this approximation. Assuming Q to 
be valid, we may reconstruct the beamformed signal using the 
Xampling schemes proposed in (121, O- 

B. Compressed Beamforming with Distorted Analog Kernels 

An obvious problem is that ^ {t;0) does not exist in the 
analog domain, and therefore may not be Xampled directly. 
We now propose a modified Xampling scheme, which allows 
extraction of its necessary low-rate samples, by sampling 
filtered versions of (pm {t) at sub-Nyquist rates. 

Since the support of ^{t;0) is contained in [0,Tb(6>)), 
where Tb {0) < T, we may define ^ {t;Oys Fourier series 
with respect to the interval [0,T). Denoting by cj the kjth 
Fourier series coefficient of ^ (t; ^), we have 

9 = ^l^ Ilo,Ts(9)) it) $ {t; e) e-^^'=^*dt, (20) 

where I[a,h) (^) is the indicator function, taking the value 1 for 
a <t < b and otherwise. Plugging the indicator function in 
(|2Q1) may seem unnecessary. However, once transforming (l2Qb 
into an operator applied directly to {(pm it serves an 

important role in zeroing intervals, which are assumed zero 
according to (O, but, in any practical implementation, contain 
noise. Substituting ^ into (l2Qb , we can write 



1 



M 
m=l 



(21) 



Fig. 4. Xampling scheme utilizing distorted exponential kernels. 



where, from Q, 



(22) 



and 



/ ^ / N / COS^ \ 

q.^mm =/[|^^|,T.(,)) W 1 + , . X 

\ (t-7^sm6') J 

f . 27r ^ 7^ - t sin 6> 1 
exp < ^^kj- ^7m ) , 

Tm{0) =rm{TB{0);e). 



(23) 



The process defined in ([2T1)-(|23]) can be translated into 
a multi-channel Xampling scheme, such as the one de- 
picted in Fig. m Each signal (pm {t) is multiplied by 
a bank of kernels {gj,m {t'',^)}f=i defined by ([23]), and 
integrated over [0,T). This results in a vector Cm = 
[ ci^rn C2,m CK,m ] • The vcctors {cmj^^i are then 

averaged in c = [ ci C2 ... ck ] , which has the desired 
improved SNR property, and provides a basis for extracting 
the 2L parameters which define ^ {t; 0). Since ^ {t; 6) satisfies 
Q, we apply a similar derivation to that outlined in SectionlTVl 
yielding 



c = 7;;HVb, 

T 



(24) 



where H is a diagonal matrix with jth diagonal element 
H{^kj), V contains e~*^^^^^ as its (j, O^h element, and 
b is the length L vector, with elements hi. The matrix V 
may be estimated by applying spectral analysis techniques, 
allowing for the vector of coefficients b to be solved by a 
least squares approach |[T9ll . Fig. [5] illustrates the shape of the 
resulting kernels gj^jji (^;^), setting 6> = and choosing two 
arbitrary values of kj . For each choice of kj we plot the kernels 
corresponding to 7 receiver elements, selected from an array 
comprising 64 elements, distanced 0.49mm apart. 



8 




t/T 



(b) 

Fig. 5. Real part of gj^m{t', = 0) for T = 210/isec and kj satisfying: (a) 
kj = 3, (b) kj=5. We assume an array comprising M = 64 elements, 
distanced 0.49mm apart, and plot 7 traces which were obtained for the 
elements indexed {mo, mo + 5, mo + 10, mo + 30}. 

VI. Simplified Xampling Mechanism 

In the previous section, we developed a Xampling approach 
to extract the Fourier series coefficients of ^ {t; 6). However, 
the complexity of the resulting analog kernels, together with 
their dependency on 9, makes hardware implementation of 
the scheme depicted in Fig. [4] complex. Here, we take an 
additional step, which allows the approximation of {cj^m}^^i^ 
and consequently {cj} -^-^^ from low-rate samples of (^rn {t), 
obtained in a much more straightforward manner. 

We begin by substituting (pm (t) of (l22l) by its Fourier series, 
calculated with respect to [0,T). Denoting the nth Fourier 
coefficient by [n] , we get: 

1 

^ ^0 ^25) 

n 

where Qj^m-e [^] are the Fourier series coefficients of 
Qj^mit^O), also defined on [0,T). Let us replace the infinite 
summation of (1251) by its finite approximation: 

N2 

^j,m = ^ (l)m[kj - n]Qj^rn-e N- (26) 

n=A/'i 

The following proposition shows that this approximation can 
be made sufficiently tight. 

Proposition 1. Assume that j'^^\iprn(f)\^ dt < oo. Then, 
for any e > 0, and for any selection (j, m; 0), there 
exist finite Ni {e^kj^m^O) and N2 {e^kj,m;0) such that 

Proof: Let I2 be the space of square- summable sequences, 
with norm ||x||2 = Let a = {0^ [kj - n]}^^_^ 

and b = < Q* ^.q [n] \ . Since Lprn (f) is of finite energy, 
a G /2. We may calculate the I2 norm of b, based on the defini- 
tion of qj^rn {t; 0) in (O, resulting in ||b||2 ^ {9) /T < 00. 



This implies that b G ^2 as well. Let bt be the truncated 
sequence b for Ni < n < N2 and zero otherwise. We may 
then write the approximation error as: 

\cj,m - Cj,m\^ = I (a, b - bt) P < ||a||^||b - btll^, (27) 

where (•, •) is the inner product defined as (x, y) = ^^XnVn- 
The last transition in (l27l) is a result of Cauchy-Schwartz 
inequality. By definition of bt and b, it is readily seen that 
l|b-bt||i = ||b||i - llbtlli. Denoting p2 = ||bt||i/||b||i, m 
becomes 

|c,,„.-c,,„P< ||a||2||b||i(l-/,2). (28) 

Since ||b||2 < 00, can approach 1 as close as we desire, 
by appropriate selection of A^i and N2. For any e > 0, there 
exists (e) < 1, such that the right side of (1281) is smaller 
than e. Selecting Ni and 7V2 for which ||bt ||^/||b||^ > (e)^ 
results in \cj^rn — Cj^m\^ < e, as required. Furthermore, setting 
an upper bound on the energy of (^rn {t), and thereby on 
||a||2, A^i and N2 may be chosen off-line, subject to the decay 
properties of the sequence {Qj,m;6> N}^_oo- * 
Using Proposition 1, we can compute Cj^m as a good approx- 
imation to Cj^rn- We now show how Cj ui 

can be obtained 

directly from the Fourier series coefficients (j)m [n] of each 

We first evaluate A^i and N2 for a certain choice of m 
and 6, such that Cj^m may be approximated to the desired 
accuracy using ([26l) . Equivalently, we obtain the minimal 
subset of (frn (^)'s Fourier series coefficients, required for the 
approximation of Cj^rn- Performing this for all 1 < j < i^, we 
obtain K such subsets. Denoting the union of these subsets 
by tvrn, we may now simultaneously compute {cj^rn}f=i ^om 
{(j)m ['^]}neKm ^ linear transformation. Define the length- 
Km vector with Ith element (j)m [h], and ki being the Ith 
element in hZrn- Using (l26l) , we may write 

= {0) (29) 

where Cm is the length-K vector with jth element Cj^m, and 
Am {0) is a K X Km matrix with elements 

^ . ^ f Qj,m-,e [kj - ki] Ni (kj) < kj -ki<N2 {kj) 
•^'^ I otherwise 

(30) 

Notice, that we have omitted the dependency of Ni and N2 on 
e, m and 6, since, unlike kj, these remain constant throughout 
the construction of A^ {0). 

The resulting Xampling scheme is depicted in Fig. [6l 
Based on [12J, we propose a simple mechanism for obtaining 
the Fourier coefficients in each individual element: a linear 
transformation, Wm, is applied to point- wise samples of the 
signal, taken at a sub-Nyquist rate, after filtering it with an 
appropriate kernel, {—t), such as the Sum of Sines. In this 
scheme, while we do need to extract larger number of samples 
at the output of each element, as Km > K, we avoid the use 
of complicated analog kernels as in Section FV-Bl Furthermore, 
as we show in Section [iXl in an actual imaging scenario good 
approximation is obtained with just a small sampling overhead. 
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Fig. 6. Xampling scheme utilizing Fourier samples of ipm (t). 



VII. Signal Reconstruction 

So far we derived our approach for extracting the param- 
eters {ti^bi}f^^ which determine ^{t;0) from sub-Nyquist 
samples, taken at the individual receiver elements. In this 
section we focus on the reconstruction of ^ {t; 6) from these 
parameters. Once ^ {t;0) is constructed for multiple values 
of 6>, a two-dimensional image may be formed, by applying 
standard post-processing techniques: first, ^{t;0)'s envelope 
is extracted using the Hilbert transform (20]; logarithmic 
compression is then applied to each envelope, resulting in a 
corresponding scanline; finally, all scanlines are interpolated 
onto a two-dimensional grid. Having obtained the parametric 
representation of ^ {t; 6), the first two steps may be calculated 
only within the support of the recovered signal. 

In Section IVII-AI we describe the reconstruction of ^ {t; 6) 
from its estimated parameters, while generalizing the model 
proposed in ©i we assume that the detected signals are 
additionally parametrized by unknown carrier phases of the 
reflected pulses, and show that the Xampling approach allows 
estimation of these unknown phases. 

In Section IVII-BI we propose an alternative approach for 
reconstructing ^ (t] 6), using CS methodology. 

A. Signal Reconstruction Assuming Unknown Carrier Phase 

Consider the signal defined in Q. Modeling a signal of 
physical nature, it is obviously real- valued, implying that ai^rn 
are real. Consequently, by ([T9l) , hi must also be real-valued. 
However, when we apply spectral analysis techniques aimed 
at solving the system formulated in (l24l) . there is generally no 
constraint that b be real- valued. Indeed, solving it for samples 
obtained using our proposed Xampling schemes, the resulting 
coefficients are complex, with what appears to be random 
phases. In fact, a similar phenomenon is observed when solv- 
ing dS]) for samples taken from the individual signals, Lprn {t), 
as proposed in [12J. Below we offer a physical interpretation 
of the random phases, by generalizing the model proposed 
in (O. The result is a closed-form solution for reconstructing 
the estimated signal, using the complex coefficients. When 
applied, a significant improvement is observed, comparing the 
envelope of the reconstructed signal, with that of the original 
signal. 

The ultrasonic pulse h (t) may be modeled by a baseband 
waveform, g{t), modulated by a carrier at frequency /q: 
h{t) = g {t) cos {ujot + /3) , where ujo = 27r/o and /3 is the 
phase of the carrier. The model proposed in (O, just like the 
one in Q, assumes the detected pulses to be exact replicas 



of h{t). However, a more accurate assumption is that each 
reflected pulse undergoes a phase shift, based upon the relative 
complex impedances involved in its reflection [21]. We thus 
propose to approximate the beamformed signal as: 

L 

^ {t; 0) = Y^ \bi I 9{t- ti) cos (a;o {t - U) + A), (31) 

Pi being an unknown phase. The jth Fourier series coefficient 
of ^ {t; 0) is now given by 



T L 

El 

1=1 



\bi\ g {t - ti) cos {uJo {t-ti) 



-^o)) < 



^=1 



-iujjti 



(32) 



where G {cu) is the CTFT of g {t) and ujj = 4^kj. 

Let g (t) be approximated as a Gaussian with variance 
and assume that kj > 0. It is readily seen that 

G {(jjj + cjo) 



We can then choose 



^o) 



-2a ojjOjQ 



5T 



SO that 



G {idj + UJo) 



< 10- 



G{ijOj - ujq) 
This allows (l32l) to be approximated as 

L 



1 



' 2T 
and additionally 

Combining ([36b and (1371) . we get 



^o) 



^H(^^k,^pie- 



(33) 



(34) 



(35) 



(36) 



(37) 



(38) 



where we define bi = |6/|e* 

Denoting by c the length K vector, with cj as its jth ele- 
ment, the last result may be brought into the exact same matrix 
form written in (l24l) . However, now we expect the solution to 
extract complex coefficients, of which phases correspond to the 
unknown phase shifts of the reflected pulses, Zbi = A — /3. 
Having obtained the complex coefficients, we may now re- 
construct ^ {t;0) according to ([311 , and then proceed with 
standard post-processing techniques. The constraint imposed 
in ([34l) is mild, considering nominal ultrasound parameters. 
Assuming, for instance, T = 210/isec, /o = 3MHz, and 
a = 630nsec, we must choose kj > 12. The requirement that 
H of ([25l) be invertible, already imposes a stronger constraint 
on kj, the jth Fourier coefficient, since H {^kj) drops below 
-3dB for I A:, - 6301 > 44. 
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B. CS Approach for Signal Reconstruction 

Throughout the previous sections, we addressed the problem 
of ultrasound signal reconstruction, within the FRI framework. 
As shown in |6|, for various FRI problems, the relationship be- 
tween the unknown signal parameters and its subset of Fourier 
series coefficients takes the form of a spectral analysis prob- 
lem. The latter is then typically solved by applying techniques 
such as annihilating filter |18 | or matrix pencil |19|. In this 
section, we consider an alternative approach for reconstructing 
the signal defined in (O, based on CS methodology H, (S). 

Assume that the time delays {ti}jL-^ in (|3T1) are quantized 
with a As quantization step, such that U = g/Ag, qi G Z. 
Using (1381) , we may write the Fourier series coefficients of 
<!> {t; 0) as: 

c,^^H(^^k,)j2bie-^'^^^^^ (39) 

Let be the ratio [T/A^J. Then ([39l ) may be expressed in 
the following matrix form: 

c ^ ^HVx = Ax, (40) 

where H is the K x K diagonal matrix with H {^kj) as its 
jth diagonal element, and x is a length N vector, whose jth 
element equals hi for j = qu and otherwise. Finally, V is a 
K X N matrix, formed by taking the set Hi of rows from an 
N X N FFT matrix. 

The formulation obtained in (l40l) , is a classic CS problem, 
where our goal is to reconstruct the A/'-dimensional vector x, 
known to be L-sparse, with L <C A^, based on its projection 
onto a subset of K orthogonal vectors, represented by the rows 
of A. This problem may be solved by various CS methods, as 
long as the sensing matrix A satisfies desired properties such 
as the Restricted Isometry Property (RIP) or coherence. 

In our case, A is formed by choosing K rows from the 
Fourier basis. Selecting these rows uniformly at random it 
may be shown that if 

K >CL{\ogN)^ , (41) 

for some positive constant C, then A obeys the RIP with large 
probability ll22ll . As readily seen from (|4T1) , the resolution of 
the grid, used for evaluating {Ujf^-^, directly effects the RIP. 
Recall that, by applying spectral analysis methods, one may 
reconstruct x from a minimal number of 2L samples, if it is 
indeed L-sparse. However, these samples must be carefully 
chosen. Using matrix pencil, for instance, the sensing vectors 
must be consecutive. Moreover, in any practical application, 
the measured data will be corrupted by noise, forcing us to use 
oversampling. In contrast, the bound proposed in (|4T1) regards 
random selection of the sensing vectors. Additionally, applying 
the CS framework, we may effectively cope with the more 
general case, of reconstructing x which is not necessarily L- 
sparse. 

VIII. Comparison between Recovery Methods 

In this section, we provide results obtained by applying 
three recovery algorithms to ultrasound signals which were 




"Z 



B 

Fig. 7. Field II simulation setup: M = 64 elements are aligned along the 
X axis with a 0.05mm kerf. The width of each element is 0.44mm. Speckle 
pattern is simulated by randomly distributing 10^ point reflectors within the 
box B. Additionally, L = 6 point reflectors are aligned along the z axis, also 
within the boundaries of the box. The pulse is transmitted along the z axis, 
and the beamformed signal is constructed along the same line. 

simulated using the Field II program (23). The evaluation 
was performed based on multiple beamformed signals, each 
calculated along the z axis (0 = 0) for a random phantom 
realization. The phantom comprised L strong reflectors, dis- 
tributed along the z axis, and multiple additional reflectors, 
distributed throughout the entire imaging medium. A mea- 
surement vector was obtained by projecting the beamformed 
signal onto a subset of its Fourier series coefficients. Finally, 
each algorithm was evaluated for its success in recovering the 
strong reflectors' positions from the vector of measurements. 
The first two algorithms which were evaluated were matrix 
pencil |[T9l and total least-squares approximation, enhanced 
by Cadzow's iterated algorithm |24|. Both algorithms may be 
considered spectral analysis techniques. The third algorithm 
was Orthogonal Matching Pursuit (OMP) (221, which is a CS 
method. 

The simulation setup is depicted in Fig. [71 We created 
an aperture comprising 64 transducer elements, with cen- 
tral frequency /o = 3.5MHz. The width of each element, 
measured along the x axis, was c//o = 0.44mm, and the 
height, measured along the y axis, was 5mm. The elements 
were arranged along the x axis, with a 0.05mm kerf. The 
transmitted pulse was simulated by exciting each element with 
two periods of a sinusoid at frequency /o, where the delays 
were adjusted such that the transmission focal point was at 
depth r = 70mm. Additionally, Hanning apodization was used 
during transmission, by applying an appropriate excitation 
power to each element. 

In each iteration, we constructed a random phantom, for 
which we simulated the beamformed signal. The phantom 
was constructed in two stages. We first created a speckle 
phantom, by drawing positions of 10^ point reflectors uni- 
formly, at random, within the three-dimensional box B = 
{(x, z) : \x\ < 25mm, \y\ < 5mm, \z — 60| < 30mm}. The 
corresponding amplitudes were also drawn randomly, with 
zero-mean and unit- variance Normal distribution. We then 
generated a signal phantom, by drawing positions of L = 6 
point reflectors, {pi}^^, with xi = yi = and zi uniformly 
distributed in the interval [35mm, 85mm). These reflectors 
were assigned identical amplitudes, which were adjusted ac- 
cording to the SNR requirement, in the following manner: for 
each of the two phantoms, we simulated the beamformed sig- 
nal, acquired along = following pulse transmission in the 
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Fig. 8. Image obtained by applying standard imaging techniques to an 
individual phantom realization. Our goal is to recover the L = 6 strong 
reflectors aligned along the z axis. 10^ point reflectors were distributed in 
the imaging plain, resulting in echoes which corrupted the detected signals. 
In the ultrasound image, these reflections are manifested in a speckle pattern. 
The phantom was calibrated such that the SNR of the beamformed signal, 
calculated along ^ = 0, defined in ( 142^ , was 15dB. 



same direction. Denoting the beamformed signal obtained for 
the first (speckle) phantom by n (t; 6> = 0) and that obtained 
for the second (signal) phantom by ^ (t; 6> = 0), we defined 



SNR = 10 logio 



Jo 



(42) 



The amplitudes of the reflectors comprising the second phan- 
tom were modified, such that (l42l) complied with the desired 
SNR value. After this calibration, we combined the two phan- 
toms into a single one, for which we generated an individual 
beamformed signal realization. The detected signals and the 
resulting beamformed signal were simulated at sampling rate 
js = lOOMHz. Since the spectrum of the detected pulses 
decayed to — 50dB at ^ 6MHz, this rate was far beyond 
Nyquist. Manning apodization was used for constructing the 
beamformed signal, by applying appropriate weights to the 
detected signals. This type of apodization may be easily 
implemented with both our Xampling schemes, by replacing 
the average in ([2T1) by a weighted one. 

Fig. [8] illustrates the method by which we simulated a 
realization of the noisy beamformed signal. This image was 
obtained by applying standard imaging techniques to an in- 
dividual phantom. We are interested in recovering the strong 
reflections aligned along the z axis. The corresponding beam- 
formed signal was corrupted by speckles, originating in the 
multiple point reflectors scattered throughout the medium. The 
phantom was calibrated such that the SNR of the beamformed 
signal along 6 = 0, defined in (l42l) . was 15dB. 

Having generated the beamformed signal, we obtained a 
measurement vector, by projecting the signal onto a subset 
of its K Fourier series coefficients, where K = 2[?7L] + 1, 
and 7^ > 1 is the desired oversampling factor. For the spectral 
analysis techniques, we chose the coefficients consecutively, 
around ko = \foT]. OMP was tested using both this se- 
lection of coefficients, and a random selection, taken such 
that H {^kj) is above — 2dB. With this selection, we obtain 
samples which are better spread in the frequency domain. 
We emphasize, that the coefficients were drawn once, for 




Fig. 9. h (t) evaluated from the beamformed signal, calculated for a 
single reflector using Field II simulator. The reflector was positioned at the 
transmission focal point. 



each choice of r]. An additional degree of freedom, using 
the OMP method, regards the density of the reconstruction 
grid, determined by N. We set N = 1860, complying with a 
sampling frequency fs = 20MHz, of order typically used in 
imaging devices. 

Recovery was evaluated based on the estimated time delays. 
These were compared to the delays associated with the known 
reflector positions, ti = 2zi/c. At the ith iteration, we 
examined, for each algorithm, all possible matches between 
actual delays {t/}^^-^, and estimated delays Of all 

possible permutations (a total number of L!), we selected the 
one for which the number of matches, achieving error smaller 
than the width of h (t), was maximal. Denoting this number by 
Sl^\ q G {1,2,3,4} corresponding to the evaluated method, 
we estimate the probability of recovery by the qth method as 



(43) 



where / is the total number of iterations, set to 500 in 
our simulation. We note that all reconstruction algorithms 
require that we first calculate H {^kj). For this purpose, we 
simulated the signal beamformed along 6> = 0, for a phantom 
which comprised a single reflector at the transmission focal 
point (x^y^z) = (0,0,70mm). We used the detected signal, 
depicted in Fig. [9l for calculating H {^kj). 

The simulation results obtained for multiple combinations 
of SNR and oversampling factor are illustrated in Fig. [Tol 
The calculated recovery probabilities are represented by gray- 
levels, where a common color-bar was used for all plots. For 
clarity, we plotted a line separating between probabilities lower 
than 0.85 and probabilities above 0.85, and a line separating 
between probabilities lower than 0.97 and probabilities above 
0.97. Of the two spectral analysis techniques, matrix pencil 
appears preferable, as it obtains high probability values over 
a wider range of SNR and oversampling. Both OMP methods 
outperformed the spectral analysis ones, with an obvious 
advantage to random OMP. 

An additional aspect which should be taken into consider- 
ation, when choosing the reconstruction method, regards the 
complexity of the Xampling hardware. Using the Xampling 
scheme proposed in [12], random selection of Fourier series 
coefficients will increase the hardware complexity: in such 
case, the sampling kernel, e.g. SoS, must be specifically 
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Random OMP 




Fig. 10. Probability of reconstruction vs. SNR and oversampling factor, 77, using four methods: (a) Total least-squares, enhanced by Cadzow's iterated 
algorithm, (b) matrix pencil, (c) OMP with consecutive Fourier series coefficients, (d) OMP with Fourier series coefficients randomly distributed, such that 



H(5 



is above — 2dB, \/kj G ^c. Signals were simulated using Field II program, where SNR is defined in (l42l . 



designed for the choice of coefficients. This is in contrast with 
the relatively simple kernel, applied for a consecutive choice 
of coefficients. On the other hand, the Xampling scheme 
proposed in 1 13] is practically invariant to the manner in which 
the coefficients are selected. 

IX. Experiments on Cardiac Ultrasound Data 

In this section, we examine results obtained by applying 
our Xampling schemes, illustrated in Figs. |4] and [6l to raw 
RF data, acquired and stored for cardiac images of a healthy 
consenting volunteer. The acquisition was performed using a 
GE breadboard ultrasonic scanner of 64 acquisition channels. 
The transducer employed was a 64-element phased array 
probe, with 2.5MHz central frequency, operating in second 
harmonic imaging mode: 3 half cycle pulses are transmitted at 
1.7MHz, resulting in a signal characterized by a rather narrow 
bandpass bandwidth, centered at 1.7MHz. The corresponding 
second harmonic signal, centered at 3.4MHz, is then acquired. 
The signal detected in each acquisition channel is amplified 
and digitized at a sampling-rate of 50MHz. Data from all 
channels were acquired along 120 beams, forming a 60° 
sector, where imaging to a depth of z = 16cm, we have 
T = 207/isec. The imaging results are illustrated in Fig. [TT] 

The first image (a) was generated using the standard tech- 
nique, applying beamforming to data first sampled at the 
Nyquist-rate, and then down-sampled, exploiting its limited 
essential bandwidth. For a single scanline, sampling at 50MHz, 
we acquire 10389 real- valued samples from each element, 
which are then down-sampled, to 1662 real-valued samples, 
used for beamforming. The resulting image is used as refer- 
ence, where our goal is to reproduce the macroscopic reflectors 
observed in this image with our Xampling schemes. 

We begin by applying the scheme illustrated in Fig. (H 
utilizing the analog kernels defined in (1231) . Modulation with 
the kernels is simulated digitally. Assuming L = 25 reflectors, 
and using two-fold oversampling, k comprises K = 100 
consecutive indices. With such selection, the corresponding 
frequency samples practically cover the essential spectrum of 
h{t). Since each sample is complex, we get an eight-fold 
reduction in sample-rate. Having estimated the Fourier series 
coefficients of <l> (t; 0), we obtain its parametric representation 
by solving (|40|) using OMP. We then reconstruct ^ {t;6) 
according to (|3T1) . that is we apply phase shifts to the modu- 



lated pulses, based on the extracted coefficients' phases. The 
resulting image (b) depicts the strong perturbations observed 
in (a). Moreover, isolated reflectors at the proximity of the 
array (z ^ 6cm) remain in focus. 

We next apply the approximated scheme, illustrated in 
Fig. [6l for every kj e k,, 1 < m < M and 0, we find 
Ni and N2 of (l26l) such that ^ 0.95. This process is 
performed numerically, off-line, based on our imaging setup. 
Consequently, we construct {Am}^^i off-line, according to 
(l3Ql) . Choosing this level of approximation, we end up with a 
seven-fold reduction in sample rate, where, for the construction 
of a single scanline, an average of 116 complex samples 
must be taken from each element. We point out that in this 
scenario, the maximal number of samples, taken from certain 
elements, reaches 133 for specific values of 6. Thus, if a 
common rate is to be used for all sensors, for all values of 0, 
we may still achieve a six-fold reduction in sample rate. As 
before, we use OMP in order to obtain ^ {t; 6)'s parametric 
representation, and reconstruct it based on our generalized FRI 
model proposed in (|3T1) . The resulting image (c) appears very 
similar to (b). 

Table IJ gathers SNR values, calculated for the beamformed 
signals estimated using both our Xampling schemes, after 
envelope detection with the Hilbert transform. The values were 
calculated with respect to the envelopes of the beamformed 
signals, obtained by standard imaging. Explicitly, let ^ {t; Oi) 
denote the beamformed signal obtained by standard beam- 
forming along the direction 6>^, i = 1,2,...,/, let ^{t;Oi) 
denote the beamformed signal reconstructed from the param- 
eters recovered by compressed beamforming along the same 
direction, and let H {•) denote the Hilbert transform. For the 
set of / = 120 scanlines, we defined the SNR as 

SNR = 10 logio 

(44) 

This calculation was repeated when reconstructing the sig- 
nals without the random phase assumption, proposed in Sec- 
tion IVII-AI For the latter case, reconstruction of a real- valued 
^{t;Oi), given complex coefficients, may be heuristically 
achieved by either ignoring the coefficients' imaginary part, or 
by taking their modulus. It may be seen that, weighting over 
all 120 beamformed signals, the random phase assumption 



EUlo\m{t;Oi))\'dt 
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achieves a relatively minor improvement (0.1-0.15dB) com- 
pared to reconstruction using the modulus of the coefficients. 
However, when examining individual signals, we observed 
that, for certain values of Oi, the improvement exceeded 1.5dB. 

TABLE I 

SNR IN [dB] of $ (t; 9) obtained with the proposed Xampling 

SCHEMES AND THREE RECONSTRUCTION METHODS 



Reconstruction Method 


Xampling Method 


Distorted Kernels 
(Fig.S 


Approximated 
Scheme (Fig.[6j 


Random Phase 


6.47 


5.89 


Real part of Residues 


4.59 


4.03 


Modulus of Residues 


6.32 


5.79 



We emphasize, that the calculated SNR values provide 
a useful measure for quantitatively comparing the different 
Xampling and reconstruction approaches. However, they are 
of smaller value when attempting to evaluate the overall 
performance of Xampling, compared to standard imaging: 
recall that our scheme is aimed at reproducing only strong 
pulses, reflected from macroscopic reflectors. The reference 
signal, on the other hand, generated by standard technique, 
already contains the additional speckle component, caused 
by multiple microscopic perturbations. A possible approach 
for evaluating the overall performance of either Xampling 
scheme, would be to examine its success rate in recovering 
strong reflections, detected by standard beamforming. For this 
purpose, we tracked the L strongest local maxima in each 
beamformed signal. If the Xampling scheme recovered a pulse 
within the range of 1.2mm from a certain maximum, we say 
that this maximum was successfully detected. Certain pulses, 
detected by Xampling, may match more than one maximum 
in the beamformed signal. In such case, we choose the one- 
to-one mapping which achieves smallest MSE. Applying this 
evaluation method to signals Xampled using our approximated 
scheme, and reconstructed with the random phase assumption, 
we conclude that the reconstruction successfully retrieves 
70.4% of the significant maxima, with standard deviation of 
the error being approximately 0.42mm. 

X. Conclusion 

In this work, we generalized the Xampling method proposed 
in (m, to a scheme applied to an array of multiple receiv- 
ing elements, allowing reconstruction of a two-dimensional 



ultrasound image. At the heart of this generalization was the 
proposal that the one-dimensional Xampling method derived 
in [12J be applied to signals obtained by beamforming. Such 
signals exhibit enhanced SNR, compared to the individual 
signals detected by the array elements. Moreover, they depict 
reflections which originate in a much narrower sector, than 
that initially radiated by the transmitted pulse. A second key 
observation, which made our approach feasible, regarded the 
integration of the beamforming process into the filtering part 
of the Xampling scheme. 

The first approach we purposed comprised multiple modula- 
tion and integration channels, utilizing analog kernels. We next 
showed that the parametric representation of the beamformed 
signal may be well approximated, from projections of the 
detected signals onto appropriate subsets of their Fourier series 
coefficients. The contribution of our schemes regards both 
the reduction in sample rate, but additionally, the resulting 
reduction in the rate of data transmission from the system 
front-end to the processing unit. In particular, our second 
approach is significant even when preliminary sampling is 
performed at the Nyquist-rate. In such a case, it allows a 
reduction in data transmission rate, by a relatively simple 
linear transformation, applied to the sampled data. 

An additional contribution of our work regards the method 
by which we reconstruct the ultrasound signal, assumed to 
obey a specific FRI structure, from a subset of its frequency 
samples. Rather than using traditional spectral analysis tech- 
niques, we formulate the relationship between the signal's 
samples to its unknown parameters as a CS problem. The 
latter may be efficiently solved using a greedy algorithm 
such as the OMR We show that, in our scenario, CS is 
generally comparable to spectral analysis methods, managing 
to achieve similar success rates with sample sets of equal 
cardinality. Moreover, working in a noisy regime, CS typically 
outperformed spectral analysis methods, provided that the fre- 
quency samples were highly spread over the essential spectrum 
of the signal. Using actual cardiac data, a relatively large 
number of reflectors was assumed. Consequently, by simply 
choosing the Fourier series coefficients consecutively, as in 
the spectral analysis techniques, we end up with the necessary 
wide distribution. However, as shown in our simulations, CS 
approach inherently allows a wide distribution of samples, 
even when the cardinality of the sample set is small, since 
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we are not obliged to unique configurations of samples. 

A final observation discussed in our work, regards the 
generalization of the signal model proposed in 1 12], allowing 
additional, unknown phase shifts, of the detected pulses. 
We show that these shifts may be estimated by appropriate 
interpretation of the extracted coefficients, without changing 
the recovery method. 

Combining the random phase assumption with our proposed 
Xampling schemes and the CS recovery method, we construct 
two-dimensional ultrasound images, which well depict strong 
perturbations in the tissue, while achieving up to seven- 
fold reduction of sample rate, compared to standard imaging 
techniques. 



Appendix A 

BEAMFORMED SIGNAL SUPPORT 

We assume h {t) to be supported on [0, A), and that the 
support of (fm {t) is contained in [0,T). The last assumption 
may be justified by the fact that the pulse is transmitted at 
t = 0, such that reflections may only be detected for t > 0. 
Additionally, the penetration depth of the transmitted pulse 
allows us to set T, such that all reflections arriving aX t > T 
are below the noise level. 

For all 1 < / < L and 1 < m < M: 



ti^m + A < T, 



(45) 



Applying the relation ti^rn = rm{ti;0), justified in Sec- 
tion IV-AI and using the fact that Tm {t; 0) is non-decreasing 
for t > we conclude that 



ti<T-' {T-A;e), 
{t; 6) being the inverse of Tm {t; 9). Explicitly: 



t—lm sin 



(46) 



(47) 



Assuming that A <C T, then, since (l46b is true for every 
1 < m < M, we may write: 



ti< min T-^ {T;0). 

l<m<M 



(48) 



This allows us to set the following upper bound on the support 
of $(t;6'): 



TbW= min r-i (T;^), 

l<m<M 



(49) 



once again, using the assumption that A <C T. From (|47]) it 
is readily seen that Tb (0) < T, since we can always find 7^ 
with sign opposite to that of sin 6, such that: 

- 72 2^2 _ 2 

^-^^■''^= T+\JZo\ ^-T^^^- '''' 

Finally, by construction of Tb {0) we see that, for all 1 < 

m<M,Tm{TB{0);0)<T. 
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